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We investigate the mixing properties of scalars stirred by spatially smooth, divergence- 
free flows and maintained by a steady source-sink distribution. We focus on the spatial 
variation of the scalar field, described by the dissipation wavenumber, kd, that we define 
as a function of the mean variance of the scalar and its gradient. We derive a set of 
upper bounds that for large Peclet number (Pc ^ 1) yield four distinct regimes for 
the scaling behaviour of kd-, one of which corresponds to the Batchelor regime. The 
transition between these regimes is controlled by the value of Pe and the ratio p = £u/£s, 
where and £s are respectively, the characteristic lengthscales of the velocity and source 
fields. A fifth regime is revealed by homogenization theory. These regimes reflect the 
balance between different processes: scalar injection, molecular diffusion, stirring and 
bulk transport from the sources to the sinks. We verify the relevance of these bounds 
by numerical simulations for a two-dimensional, chaotically mixing example flow and 
discuss their relation to previous bounds. Finally, we note some implications for three 
dimensional turbulent flows. 
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1. Introduction 

Mixing of scalar fields is a problem that is crucial to several environmental issues as well 
as engineering applications. In many situations the underlying flow is spatially smooth 
and divergence-free while mo lecular diffu sion is usually much weaker than the stirring 
strength of the flow (see e.g. Are! ( 20021 )). Notwithstanding the apparent simplicity of 



the flow, its effect on the scalar field can be rather complex: A simple time-dependence 
is often sufficient for the flow to be chaotically mixing in whi c h case the gr adients of 
the sc alar flclds are greatly amplifled (jArej (|l984[) : lOttinol (|l989t ):[Ottl (|l993l )). lBa"tchelor 



( 1959( ) recognized that this amplification is responsible for the rapid dissipation of any 



initial scalar inhomogeneity and thus the efficiency at which a scalar is mixed. 

In the continual presence of sources and sinks, a statistical equilibrium is attained 
in which the rate of injection of scalar variance balances the rate of its dissipation. In 
this case, the most basic way to measure the flow's mixing efficiency is to consider the 
equilibrium variance of the scalar: the lower its value, the better mixed is the scalar 
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field. iThiffeault etai\ ([200^. derived a rigorous low er bound for the scalar variance that 



was further e nhanced bv IPlasting fc Young ( 20061) using the scala r dissipation rate as 



a constraint. Doering fc Thiffeaultl ( 20061 ) and Shaw et al. ( 2007 ) derived bounds for 



the small- and large-scale scalar variance (respectively measured by the variance of the 
gradient (|V6'p) and the anti-gradient (|V~^6'p) of the scalar field 9, where (•) denotes 
a space-time average defined in equation (|2.2p ). This set of bounds have successfully 
captured some of the key parameters in the fiow and source-sink distribution that control 
the scalar variances. Their general applicability means that they can be used to test 
theoretical predictions of scalar mixing for various flow and source-sink conflgurations. 
This is especially useful for high-Peclet flows (Pe ^1) for which analytical solutions are 
difficult to obtain while high-resolution numerical simulations can become prohibitively 
expensive. However, the bounds on the variance of the scalar and its gradient do not 
depend on the gradients of the velocity field and in many cases, can be realized by 
uniform flows. They therefore do not capture the effect of stirringjj. These bounds are 
then relevant when the mixing of a scalar is mainly controlled by the process of transport 
from the sources to the sinks. 

Motivated by the apparent lack of control of the stirring process, we here focus on the 
characteristic lengthscale, ^rf, at which the scalar variance is dissipated, or cquivalcntly its 
inverse, the dissipation wavenumber, = Its value, should, within a suitable range 
of paramete rs, be directly rel ated to the Batchelor lengthscale, . The latter lengthscale, 
obtained in Batchelor (Il959l) . describes the effect of stirring on the spatial structure of 



the scalar field. 

We here examine the behaviour of k^ for different values of the control parameters, 
Pe and p, where p denotes the ratio of the characteristic lengthscale of the velocity, £u, 
and that of the source field, ig. After formulating the problem in section [51 we next seek 
a set of upper bounds for fc^ (section [3]). In section |4l we investigate the behaviour of 
these bounds as p varies. We find that, in the high-Peclet limit, the behaviour of kd is 
characterized by four distinct regimes, one of which corresponds to the Batchelor regime. 
The use of homogenization theory implies a fifth regime for kd- In section [5l we examine 
the relevance of the bounds by performing a set of numerical simulations for a renewing 
type of flow. We conclude in section [6l 



2. Problem formulation 

The temporal and spatial evolution of the concentration, 9{x,t)^ of a passive scalar, 
continually replenished by a source-sink distribution, is given by the forced advection- 
diffusion equation. Its general form, expressed in terms of dimensional variables, is given 

by 

dte{x, t) + u{x/tu, t) ■ \/9{x, t) = kA9{x, t) + s{xm, (2.1) 

where k is the molecular diffusivity, u{x/iu,i) is an incompressible velocity fleld (i.e. 
V • M = 0) and s{x/£s) is a steady source field. Both u{x/iu,t) and s{x/£s) are spatially 
smooth (i.e. |Vs|, |(Vit)ij| < oo), respectively varying over a characteristic lengthscale £u 
and £s that can be taken to be the smallest (persistent) lengthscale in the corresponding 
fields. They are prescribed within a domain, n, that we take to be a d-dimensional box 
of size L on which wc apply either periodic or no-flux boundary conditions. This way, the 
boundaries can not generate any additional variability in the scalar field. The amplitude 

f The dependence on the velocity gradients only appears in the lower bound for the l arge-scale 
varian ce (|Shaw et al\ (|2007l )) and its decay rate in the case of no sources and sinks (|Lin et al\ 

dm!)). 
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of the velocity and source field is respectively measured hy U = {u ■ u) and S 
where (•) represents a space-time average such that 



1 



(/)= lim — - / dt dx f{x,t), (2.2) 

and Vh denotes the volume of the domain. Without loss of generality, we can assume 
that the spatial averages of 6{x,Q) and s{x) are both zero (where negative values of s 
correspond to sinks for 9) so that 9{x,t) eventually attains a statistical equilibrium with 
(9) = 0. 

We are here interested in the processes that control the mixing of 9{x,t) and how 
these depend on two non-dimensional parameters associated with equation (j2.ip . The 
first parameter is the Peclet number, Pe, defined as 

Pc = Ueu/n, (2.3a) 

which describes the strength of stirring relative to molecular diffusion. The second pa- 
rameter is the ratio, p, of the velocity lengthscale, £u, to the source lengthscale, is, defined 
as 

P = in/is- (2.36) 

There are many ways to quantify mixing. The simplest perhaps measure is given by 
the long-time spatial average of the scalar variance, which for (9) = 0, reads 

= (9^). (2.4) 

A scalar field is well-mixed when its distribution is nearly homogeneous i.e. has a value of 
a that is small. Conversely, a badly-mixed scalar distribution is one that is inhomogencous 
i.e. has a large value of a. 

The large-scale scalar variance introduced by the source at is is transferred into small- 
scales where it is dissipated by molecular diffusion. This transfer is greatly enhanced by 
the amplification of the scalar gradients induced by a stirring flow. The average rate at 
which the scalar variance is dissipated is given by 2% where: 

X^K{\Vef). (2.5) 

Wc can now define the dissipation lengthscale, id, as the average lengthscale at which the 
scalar variance is dissipated. Let the dissipation wavenumber, k^, denote the inverse of 
id- Then, id and kd are given by 



By construction, the dissipation scales (|2.6|) characterize the spatial variation of the scalar 
field and as such, provide an alternative way to quantify mixing. 

The dissipation wavenumber is related (although it is not always equ al) to the d iffusive 
cut-off scale of the 0-spectrum. For a freely decaying scalar (i.e. s = 0). iBatchelo r (1952i) 
estimated this cut-off lengthscale to be independent of the initial configuration of the 
scalar field with 

u - Vp~c ^ ' 

where i^ stands for Batchelor's lengthscale. Being independent of the source properties, 
£g can be used as a reference to which the value of kd can be compared for varying values 
of p and Pe. 



4 



A. Alexakis and A. Tzella 



Multiplying equation (|2.1|) by 9 and taking the space-time average ()2.2|) gives the 
following integral constraint for 

X = {Os). (2.8) 

Thus, the average rate at which scalar variance is injected by the source at is equal 
to the average rate at which scalar variance is dissipated by molecular diffusion at small 
scales. Using the integral constraint (|2.8p . it is then straightforward to show that fc^ and 
a are intimately related. In particular, 

les) cr2 

a = - — - X — 

X 

= i0,s -, where 6,, = (2.9) 

^6(,s expresses the correlation between the scalar and source fields and takes the values 
between ^ ^6»,s =^ 1- For fixed values of 5* and k, there exist two ways to reduce the 
variance of the scalar: The first one relies on minimizing the correlation ^e_s while the 
second one relies on maximizing the value of kd- Minimizing the correlation ^g^g can be 
achieved by choosing a flow that rapidly transports fluid parcels from a source region 
(s > 0) to a sink (s < 0). In this configuration, the flow is not necessarily a stirring 
flow; a uniform flow can be just as efficient in reducing ^e,s (see Thiffeault fc Pavliotid 
( 20081 ) where the importance of efficient scalar transport from the sources to the sinks 
is highlighted for optimal mixing). The flow process that suppresses the scalar variance 
is in this case the process of transport. The second way to reduce the scalar variance 
is to increase the value of kd- This increase can be achieved by choosing a flow that 
rapidly stretches fluid parcels so that the magnitude of the scalar gradients are greatly 
amplified. The flow process that suppresses the scalar variance is in this case the process 
of stirring. Thus, information about either ^g^s or kd can provide us with some insight 
on the mechanisms involved in the reduction of a. 

In the next session we focus on bounding the value of kd- 



3. Upper bounds for the dissipation wavenumber 

3.1. Previously derived results 

Proper manipulation of the forced advcction-diffusion equation p.l[) leads to a number 
of constraints that can be employed to deduce a set of upper and lower bounds for the 
mixing measure s under consider a tion. A first integral constraint is given by equation 
(|2.8p . Following Thiffeault et al\ (2004), a second integral constraint can be obtained 
by multiplying equation by an arbitrary, spatially smooth 'test field', i^ix), that 

satisfies the same boundary conditions as 6{x). Space-time averaging and integrating by 
parts leads to 

(Ou-Vi:) + K{eA^) ^ -{si:). (3.1) 

Choosing ip — s we first apply the Cauchy- Schwartz inequality on equation p.ip to 
isolate a. We then use Holder's inequality which leads to the following lower bound for 
the variance a: 

CT^ r, (3-2a) 

C/sup,|Vs| +k(|As|2)4' 

Cl + Pc pC2 
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where ci and C2 are non-dimensional numbers that only depend on the 'shape' of the 
source field and not on its amplitude or characteristic lengthscale. Explicitly they are 
given by 

ci and C2 = , (3.2c) 

where the hat symbol signifies differentiation with respect to x = x/ig- Note that for 
Ci and C2 to remain 0(1), is needs to represent the smallest length of variation in the 
source field. Using expressions (|2.9p and (|3.2r b')'). we obtain the following upper bound 
for kl: 

kj^ — ^^{c^+Pc-'pc2). (3.3) 

Thus, for sufficiently large Pcclet number, the upper bound for kd is determined by the 
magnitude of U/£s, the typical timescale associated with bulk scalar transport from the 
sources to the sinks, relative to n, the molecular difFusivity. Once normalized by the 
Batchelor lengthscale p.7p . expression p.3p becomes: 

kjil^p{ci+Pc-'pc2). (3.4) 



Both bounds (|33) and were first derived in iThiffeault et ali (|2004D . 

3.2. A new upper bound 

A new upper bound for fc^ can be obtained by considering the spatial and temporal 
evolution of the gradient of 9, 

dtVe + u ■ V(V6i) = kAV9 - (Vtt)^V6i + Vs, (3.5) 

where the upper index T stands for transpose and [(V^t)^V^^],; = X]j=i(^i'^i)^i^- The 
average rate at which the variance of the scalar gradient is dissipated is 2ri where ry is 
defined by 

rj^K{\Ae\'). (3.6) 

Multiplying equation (|3.5p by V6' and taking the space-time average (|2.2p gives the 
following integral constraint for 77: 

77 = -(v6'(VM)''"°ve') + (ve'-Vs), (3.7) 

where the tensor (Vm)^J™ = ^[{\/u)ij + {Vu)ji] is the symmetric part of the velocity 
gradient tensoiU. Using Holder's inequality, the first term in (|3.7p is bounded by: 

|(Ve'(Vtt)^y'"V6l)| sC sup |n(Vit)^y'"n|(|V6l|2), (3.8a) 

where n is a unit vector so that \n\ = I. Integrating by parts the second term in (j3.7p 
and using the Cauchy-Schwartz inequality results in: 

|(V0- Vs)K (T(|As|2)i (3.8&) 

Combining the two bounds in p.8p leads to the following upper bound for the dissipation 
rate of the variance of the scalar gradient: 

U X 

t Note multiplication of equation H3.5P by a 'test field', ij^ix)^ yields constraint (|3.ip . 
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where C2 and C3 are non-dimensional numbers that depend on the shapes of the source 
and velocity field, respectively. C2 was previously defined in equation p.2c[) and C3 is 
defined by 

C3 = ^sup|n(Vur'°n|, (3.10) 

where the tilde symbol signifies derivation with respect to x = x/£u- Note that for C3 
to remain 0(1), iu needs to represent the smallest persistent length of variation in the 
velocity field. 

The upper bound for 77 in equation p.9p can serve to bound kd by observing the 
following inequality that relates x^ and 77: 

K\{0A9)\i^<7y/l^, (3.11) 

obtained by partial integration and application of the Cauchy-Schwartz inequality on the 
definition of x in equation p.Sp . Using the definition (|2.6p of kd and the square of (|3.1ip 
we then have 

kj ^ XP- {3A2a) 

U rv 



where the bounds (|3.9p on rj and p.2p on a were employed in order to deduce the last 
inequality. The above quadratic inequality in k^ yields the following upper bound for fc^: 



kjel < ic3 + iy^ci+VPc-i(ciC2 + pc2pe-i), (3.13) 

where as before, k^ is normalized by the Batchelor lengthscalc (12.71) . 

Bound (|3.9p can further be improved for the particular case of a monochromatic source 
i.e., a source that satisfies the Helmholtz equation: 

As = -cafc^s. (3.14) 

It follows that |(V0- Vs)| = C2k1{9s) = cafc^Xi where the latter is directly obtained using 
the integral constraint (|2.8p . Substituting in equation p.7p . bound p.9p becomes 

f^s^ + ca^'jx. (3.15) 



From constraint (|2.8p . % ^ US' and thus equation p.lSp provides a better bound for 77 
than equation p.9p . Using this inequality, equation p. lip leads to 

klll^c^ + C2p'Ve-\ (3.16) 



4. Different regimes 

Figure [T] shows the behaviour of the two bounds, given by Eqs. p.4p and (|3.13p . for 
various Peclet numbers, as a function of p. For small Peclet number (Pc < 1), bound 
p.l3p does not improve bound (|3.4p since for all values of p it is either greater or similar 
to bound p.4p . However, as the Peclet number increases beyond 0(1) values, the process 
of stirring becomes increasingly important and expression (|3.13p can significantly improve 
the upper bound for A:^^^ . This improvement depends on the value of p. It is only for 
values of p ^ 0(1) that bound p.l3p becomes smaller than bound p.4p and thus a better 
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iD'^ r 

10"^ 10"' 10° 10' 10^ 10^ 

p 

Figure 1. The upper bounds (|3.4|l (dashed hue) and (|3.13p (sohd hne) plotted as a function 
of p for five different values of the Peclet number: Pe = fO^'^, fO'\ fO^, fO®, and in the limit 
of Pe — >■ oo (the constants Ci, C2 and cg are given in equation ([S3])). For Pe fO'', the upper 
bound p.4p remains nearly invariant within the plotted domain. 



upper bound for Thus, in the high-Pcclet limit (Pe ^1), the two bounds capture 

different regimes of mixing that we now describe. 

We first focus on p ^ 0(1). The three terms inside the square root in equation (|3.13p 
give rise to three different power-law regimes for the behaviour of the upper bound of 



4.1. Regime I 

For p » Pe, the last term inside the square root in equation p.l3p dominates. Thus, 
k'^i'i ^ C2/5^Pe~"^ whence. 



kd^^, for Pe«p, (4.1) 

where sub-dominant terms have been dropped. In the case of a monochromatic source, 
the validity of this regime extends to p ^ \/Pc. 

For this range of values of p, the flow is nearly uniform with respect to the source while 
diffusion acts faster than transport. As a result, the scalar variance that is injected by 
the source is directly balanced by diffusion. Thus, to first order, the effect of the flow can 
be ignored from where we obtain that « with c another non-dimensional number 
defined as c = (|V~^sp)^/S'. Note that for a monochromatic source, x = KC2k^a^ and 
thus bound (14.11) is saturated. 



4.2. Regime II 

For Pe3 <c p <C Pe, it is the second term inside the square root in equation p.l3|) that 
dominates. In this case, k'^i'^^ ^ y/ciC2P^/Pe and thus the following applies for k^: 

I / U( \ ^ 1 
fed ;^ ( ciC2-^ j , for Pc3</,<Pe (4.2) 

where sub-dominant terms have been dropped. The flow continues to be slowly varying 
for these values of p. However in this case, diffusion is not the only dominant process: 
the time of transport between the sources and the sinks becomes important. Bound 
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(|4.2p reflects this importance in its dependence on Pe/p, the ratio of times of diffusion 
and transport between the sources and the sinks. At the same time, the non-trivial Pe- 
dependent scaling of bound (|4.2p can not be deduced by the balance of only two processes 
(as is the case for Regime I). This scaling is likely to be related to the formation of 
boundary layers within which the scalar variance is large. Their generation is associated 
with regions in which the continual inj ection of scalar va riance cannot be suppressed by 



sweeping across the sources and sinks. IShaw et al\ (|2007l ) examined the case of a steady, 



uni-directional shear flow and a monochromatic source from where they obtained that for 

1 2 

Pe ^ 1, kd ^ Pc^ p~ 3 kg. Nevertheless, we here find that Regime II is absent in the case 
of a monochromatic source. Whether the scaling suggeste d by bound (14.21) i s realized by 
more complex flows and source functions than the one in IShaw et al\ (j2007[ ) or if bound 
(|3.13p can be improved remains an open question. 

4.3. Regime III 

The third regime appears for 0(1) ^ p <^ Pe^. In this case, the first term inside the 
square root in equation p.l3p dominates and bound p.l3p becomes ^ C3. Thus, in 
this regime, the bound for fc^£^ implies that kd and are ii iversely pr o portio nal to each 



tnis regime, the bound tor fc^t^ implies that Kd ana tg are i nversely pr o portio nal to each 
other. This relation corresponds to the prediction made in iBatcheloij (|l959f ). It follows 



kd^\lc3-^, for 0(l)s$p«Pe*, (4.3) 



that 



where sub-dominant terms have been dropped. Note the dependence of equation (|4.3p 
on the stirring timescale, iu/U. It is therefore clear that in this regime, the dissipation 
wavenumber is governed by the balance between the processes of diffusion and stirring. 
Note that for a monochromatic source, this regime appears for 0(1) ^ p <C Pe^ . 

4.4. Regime IV 

When p ^ 0(1), the characteristic lengthscale of the source becomes larger than that of 
the velocity field and bound p.4p becomes relevant. In this case, < cip and thus 



kd^xlc^^, for p^O{l). (4.4) 

Thus, in this regime, both the processes of transport between the sources and sinks and 
diffusion control the behaviour of the dissipation wavenumber. 

4.5. Regime V 

Although not captured by the two bounds, a fifth regime is expected to appear when the 
characteristic lengthscale of the fiow is much smaller than that of the source (p <C 1). In 
this case, the large-scale solution to equation (j2.1|) is well- approximated by d{x,t) that 
satisfies the following equation: 

dt9 = W-K-Wd + s, (4.5) 



where an effective diffusion operator has replaced the advective term in equation (|2.1|) . 
The effective diffusivity tensor, K, can be written as 

K^nil + Kr), (4.6) 

where I is the identity tensor and Kt is a (non-dimensional) tensor that represents the 
enhancement of the diffusivity due to the flow. It thus follows that for this range of values 
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of p, the dissipation wavenumber can be approximated by 



X 



/3< 1. 



(4.7) 



This approximation is obtained using « X ~ (•s^) and muhiplying equation (j4.5p 
by 9 and space-time averaging to estimate {s9) . 

The coefficients of Kt can be rigorously obtained within the framework of homog- 
enization theory in which multi-scale asymptotic methods are employed in order to 
de rive the large-scale effect of the s mall-scale velocity f i eld (fo r derivation see review 
by iMajda fc Kramer ( 1999f ) and also Kramer fc Keating ( 2009f ) in which the case of a 
continuously replenished scalar is examined). In general, the coefficients of Kt depend 
on the value of Pe with ||-K^t|| ^ Pe", where the exponent a depends on the type of 
flow under consideration. For shear flows (Taylor transport), a = 2; for globally mixing 
chaotic advection flows, a = 1; for cellular flows with closed field lines, a = 1/2 (see 
Majda fc Kramen (|l999l )). Thus, depending on the value of a. 



whence, 



.^Pe"- 



(4.8) 



(4.9) 



For fixed value of Pe, the above scaling increases faster in p than the bound for k-^l^ in 
Regime IV. It follows that the validity of the asymptotic result (|4.9|) is constrained by the 
upper bound p.4p . For sufficiently high Pcclet values, this is the case when p < 0(Pe^~"). 
Based on this argument, the scalings (|4.8|) and ()4.9p are expected to be valid at most 
when 

p<min{l,Pe^-"}. (4.10) 

In general, the range of vali dity of the homogenization theory is limited to p ^ Pe~^ (see 
Kramer fc Kea ting (2009); iLin et all tOldi). The re l evanc e of Pe"^ was shown to be 
true for the mixing measures of iDoering fc Thiffeault (120061 ) . calculated for a particular 



class of steady flows (wit h a = 2) in Lin et al. ( 2010l ) and for a family of steady flows 
of various values for a in iKeating fc Krameij (|201of r For chao tic flows however (a = 1) 
the predictions of homogenization theory have been shown in Plasting fc Young ( 20Q6[ ) 
to be surprisingly accurate even for p ~ 0(1)- 



5. Numerical simulations for a representative flow and source 

We now examine how close the bounds are to the dissipation wavenumber, obtained 
from the solution of the forced advection-diffusion equation (|2.ip . To that end, we per- 
form a set of numerical simulations for a passive scalar, advecte d by a renewing chaoti c 
advectio n flow, the w idely employed alternating sine flow (e.g. iPierrehumbertI (jl994l ): 
lAntonsen et all (|l996} )). This flow is explicitly given by 



U{x/iu,t) = 



Q{T/2-t mod t) %/2t/sin(y/^^ 
e(t mod t-t/2) V2UsuY{x/in 



h) 
b2) 



(5.1) 



where Q{t) is the Heaviside step function defined to be unity for t ^ Q and zero otherwise. 
(pi and 02 are independent random angles, uniformly distributed in [Q,27r], whose value 
changes at each time-interval r in order to eliminate the presence of transport barriers 
in the flow. This way the flow is globally mixing. The alternating sine flow is isotropic 
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and homogeneous in the sense that 

M = ^U%,. (5.2) 

For this flow, the Strouhal number St can be defined in terms of the stirring timescale 
and the correlation timescale, r: 

Si.^. (5^3) 

We choose a monochromatic source field that is given by 

s{x/£,) = 2Ssm{x/es) sin(?//£s). (5.4) 



This source field satis fies equation (13.141) and thus the two relevant bounds are (|3.4 



and (|3.16p . Note that Plasting fc Young ( 2006[ ) showed that for this particular set-up 



the choice of V' = s in constraint (|2.8p is an optimal one for the variance. We take the 
domain to be a doubly periodic square box whose size L is equal to the largest of the two 
spatial lengthscales L = 27r max(^„, £s). For this flow and source fields, the coefficients 
ci, C2 and C3, defined in Eqs. <\3.2c^ and p.lOp . are given by: 

ci = 2V2, C2 = 2, C3 = V2. (5.5) 
In the high-Peclet limit (Pe ^ 1), the effective diffusivity tensor in equ ation (14.61) ca n 



be calculated from the s ingle-particle diffusivity of the velocity field (see iTavloij (|1921I ) 



Maida fc Kraniei] (1999)). For flow (|5.1l) . the enhancement diffusivity tensor, Kt, is 



found (see also Plasting fc Young] (|20G6f )) to satisfy 



C/^T TrPe 
Kt = —I = ^^J. (5.6) 

Employing equation (|4.9p . we can derive the following prediction for the dissipation 
wavenumber: 



We solve the forced advection-diffusion equation \2.1\ for flow (|5.ip and source (|5.4p 
using a pseudo-spectral method with resolution of up to TV = 8192 grid points in each 
direction. We consider different values of the two control parameters, p and Pe. We first 
focus on two values of Pe: Pe = 3.5 x 10^ and Pe = 1.4 x 10^ and keep St fixed with 
St = 1. For the first value of Pe, p varies in powers of 2 between 1/32 and 128. The second 
value concerns larger values of p, varying between 1 and 128. In all simulations, the grid 
size is chosen to be smaller than the Batchelor lengthscale, 1^. Thus, N > L/l^. We let 
the simulation evolve in time until a well-observed, statistically steady state is reached. 
The time-averages of all quantities of interest are thereafter calculated over several time 
periods r. 

5.1. Scaling regimes 

Figure [2] compares the two theoretical upper bounds, (|3.4p and (|3.16p . with the numerical 
values for /c^^^. Also shown is the prediction for fc^^^, obtained from homogenization 
theory. The two upper bounds combined with the prediction of homogenization theory 
capture the non-trivial dependence of on p. In particular, the theoretical curves and 
the numerical results share, for similar range of values of p, similar slopes. 

However, the different scaling regimes associated with the bounds are more difficult to 
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Figure 2. Numerically obtained values for fcl^^ plotted as a function of p for Pe = 3.5 x 10^ 
(diamonds) and Pe = 1.4 x 10^ (squares). The values are obtained from a set of simulations 
for flow (IS.lfl and source (|5.4p for which St — 1. Also plotted for comparison the upper bounds 
Ip^it (dashed line) and (pJ6)) for Pe = 3.5 x 10^ (sohd line) and Pe = 1.4 x 10^ (dash-dot line). 
The dotted line on the left shows the prediction (|5.7p of homogenization theory while the two 
dotted lines on the right show the diffusive scaling associated with Regime I. The filled symbols 
mark the simulations associated with figures |4l[a-c). 



discern. This is not surprising since for each power-law to clearly appear, p needs to vary 
by at least an order of magnitude. This is numerically prohibitive, especially for p ^ 1 
in which case N > -x/Pe/p. At the same time, for the chosen flow (|5.ip and source (|5.4p . 
it is not clear that the dissipation wavenumber will, in each of the regimes, scale like the 
bound. 

Still, in figure [2] we see that for p <C 1, the p^-dependent prediction (|5.7p of the ho- 
mogenization theory i.e.. Regime V, is in good agreement with the numerical results. 
As p increases to 0(1) values. Regime III becomes relevant: the slope of fc^^^ decreases 
significantly with fc^^^ becoming nearly constant. This is particularly true for the sim- 
ulations corresponding to Pe = 1.4 x 10^ for which Regime III extends to a larger range 
of values of p. For simulations with Pc = 3.5 x 10'^, Regime III is limited to a smaller 
range of values of p and a transition to the diffusive Regime I appears, as demonstrated 
in the figure. Note that, as expected, the bound in Regime I is saturated (see discussion 
in section |4T|) . 

Although the homogenization prediction (|5.7|) provides a good description of kd at 
small p, its dependence on the Strouhal number suggests that the range of validity of 
Regime V can be limited. An estimate for the validity range of Regime V can be obtained 
by calculating the point of intersection between (|5.7p and p.4p . Thus, for flow (|5.ip and 
source (|5.4p . equation (|4.10p becomes 

g 

/?< min{l,-S't}. (5.8) 

TT 

According to equation (|5.8p , we expect that as the value of the Strouhal number decreases 
below 0(1) values, the transition to Regime V will occur at increasingly small values 
of p. This expectation is reflected in the numerical values for that are shown in 
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Figure 3. Numerically obtained values for fc^^g for various p and three values of St: St — 1 
(diamonds), St — 1/4 (triangles), St — 1/16 (squares). Values obtained from a set of simulations 
for flow (|5.ip and source (|5.4p for Pe = 3.5 x 10^ . The dashed and solid lines respectively indicate 
the upper bounds (|3.4I) and H3.13[l while the dotted lines indicate prediction 1)5. 7[) obtained from 
homogenization theory. 
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figure [31 obtained from a set of simulations for St ~ 1, St = 1/4 and St = 1/16 and 
p ^ 1. Note how closely prediction (|5.7p matches the numerics. At the same time, for 
min{l, ^St} ^ p ^ 0{1), the numerical results all collapse to the same power-law regime 
whose exponent /3 ~ 0.5 lies in between the corresponding exponents associated with 
Regimes III and IV. We anticipate that further decrease in the values of St and p, will 
bring out the p-dependent scaling of Regime IV. 

5.2. Spatial structures 

It is interesting to relate the variation of the dissipation wavenumber with the different 
spatial structures of the scalar field that result from varying p. Figures SJa-c) show 
three snapshots of 9{x,t) obtained for three different values of p, chosen for their clear 
representation of the different spatial structures that can be obtained. Also shown in the 
same figures are the time-averaged variance spectra, af. , defined in terms of 9q {t) , where 

eix,t)^j:Jg{t)e''-^: 



1 



E K(.t)\'dt, (5.9) 



To — Xi It 

'"'^^ fcs;|q|<fc+l 

where Ti is chosen to be sufficiently large for a steady state to have for long been estab- 
lished and Ta - Ti > r. 

FigurelUa) displays the case of p = 1/16 and Pe ~ 3.5 x 10^. This set of values corre- 
sponds more closely to the homogenization Regime V (see also figure [2]). The scalar field 
(left panel) is essentially a superposition of a large-scale component that is proportional 
to the source field and a small-scale component that is generated by the stirring veloc- 
ity field. This superposition is clearly depicted in the spectrum of the variance (right 
panel) in which we observe that the majority of the scalar variance is concentrated on 
a single wavenumber which corresponds to the characteristic wavenumber of the source, 




(c) p = 32 and Pe = 1.4 x lO'^ 



Figure 4. Snapshots of the scalar fields (left panels) and their normalized variance spectrum 
(right panel) for various values of p and Pe. The dotted line shows Batchelor's power-law 
prediction for the spectrum. The vertical dash-dot line shows the location of the Batchelor 
wavenumber = 1/lg . The vertical dash-dot-dot-dot line shows the location of the dissipation 
waver number, kd- For figure[4l^c), the dashed line in the variance spectrum (right panel) shows 
the fc"*"^ scaling that is related to the 5-like structures shown in the left panel. 
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ks = l/£s- Small-scale fluctuations are present for fc^ < k < ku, where fc„ = l/^„ is the 
characteristic wavenuniber of the velocity and k^ = 1 is the Batchelor wavenumber 
beyond which the spectrum falls off exponentially. The spectrum associated with the 
small-scale fluctuations exhibits a power-law scaling with an exponent that is somewhat 
smaller than the Batchelor prediction of —1. The difference between kd and fc^ arises be- 
cause of the large concentration of scalar variance at small wavenumbers which shifts the 
value of fcrf to values smaller than fc^ . With decreasing p, the amplitude of the large-scale 
variance increases and so does the difference between kd and k^ . 

FigureHl^b) displays the case of p = 1 and Pe = 1.4 x 10^ that most closely corresponds 
to Regime III. For this set of values, we observe the classical filamental structures that are 
obtained when stirring dominates. In this case, the field has no memory of the functional 
form of the source. For fc < fc^ , the variance spectrum is characterized by a clean power- 
law that behaves in agreement with the Batchelor prediction. 

Figure IH^c) displays the case of p = 32 and Pe = 1.4 x 10^ which corresponds to 
the region of transition between Regime I and III (recall that for source (|5.4p there is 
no Regime II). In this case, the flow is slowly varying with respect to the source. As 
a result, the scalar variance is, in the bulk of the domain, controlled by the process 
of sweeping between the sources and the sinks and thus its value is small. The only 
exception is a small number of isolated thin, boundary layers within which the scalar 
field is highly varying. These thin layers are formed in regions where the background flow 
is nearly stagnant so that the continual injection of scalar variance cannot be suppressed 
by the process of sweeping. For our particular flow ()5.ip the regions of zero velocity are 
lines that, depending on time, are either vertical or horizontal. Thus, in each half of 
the period (horizontal/vertical) thin la yers of alternating sign of 6 are formed. These 



layers are similar to those obtained in IShaw et al\ (j2007l ) for a steady, sine flow. The 



(horizontal/ vertical) thin layers that are formed in the first half of the period are then 
stretched in the second half at the same time as new (vertical/horizontal) thin layers are 
formed. The formation of these highly-varying, thin layers yields two power-law scalings 
for the variance spectrum: For small wavenumbers (fc < fc^), the spectrum has a positive, 
power-law behaviour, given by cr^. ~ k. This behaviour can be deduced by noting that, for 
scales much larger than the source lengthscale, these structures are "(5-like" . Conversely, 
for large wavenumbers - albeit smaller than the Batchelor wavenumber (kg < k < k^), 
the Batchelor spectrum is recovered. 



6. Conclusion 



In this work we obtained a set of upper bounds for the dissipation wavenumber. kd, of 
a continuously forced scalar field that is stirred by a spatially smooth velocity field. We 
focused on the dissipation wavenumber because it provides a measure for the enhance- 
ment of mixing due to the process of stirring. Unlike the freely decaying case in which 
stirring is the only mechanism for efhcient mixing, in the forced case, transport can be 
as effective as stirring. This is clear from equation p.9p in which it is easy to see that 
the scalar variance can be reduced either by increasing the dissipation wavenumber or 
by decreasing the correlation between the scalar and source fields. 

Previous investigation s have considered a number of mixing m easures for which a set of 



bound s were derived (seelThiffeault et a/.l (|2004l ): lPlasting fc Youn g ( 20061 ) : lDoering fc Thiffeauh 



(|2Q0d ): IShaw et al\ (120071) 1. ^lowever, these bounds do not always distinguish between 
the processes of stirring and transport. In particular, the bounds on the average variance 
of the scalar (0^) and its gradient {\\/9\'^) do not explicitly depend on the velocity gradi- 



Bounding the scalar dissipation scale for mixing flows 



15 



Regime 
I 

II 

III 

IV 



V 



s; ciU/k£, 
~ p^Pe"-V4 



range of validity 

p > Pe 
p > Pei 

Pes < p < Pe 
(absent) 

<P<Pe* 
<P<Pe^ 

minjl.Pe^-"} <p<0(l) 

p < min{l,Pe^""} 



type 

difFusion-dominated regime 
(monochromatic source) 



(monochromatic source) 

Batchelor regime 
(monochromatic source) 



homogenization regime 



Table 1. The four regimes deduced from bounds p.4p and (|3.13|) and from homogenization 
theory. Al so noted the case of a monochromatic source satisfying equation (|3.14l) (adapted from 
review bv lThiffeaulj ((MH )). 



ents. Thus, the effect of stirri ng is not captured . As a result, the bound for kd docs not 
follow the scaling predicted in iBatchelor (I1959I) . 

With the aid of an additional constraint, we here derived a new upper bound for 
which, within a range of values of Pe and p, is, up to a constant, equal to the inverse 
of the Batchelor lengthscale, £~^. The process of stirring is thus reflected in this bound. 
For large Peclet values, both the previous and the new bound become important, with 
the new bound significantly improving the previous bound for p > 0(1). The scalings 
associated with these bounds suggest four different regimes for kd- The use of homoge- 
nization theory implies a fifth regime. The most interesting, perhaps, behaviour occurs 
for p ~ 0(1). For these range of values of p, the scaling suggested by the upper bounds 
for kd transitions from a behaviour controlled by transport to a behaviour controlled by 
stirring. A summary of our results is provided in table [TJ 

We tested the relevance of our theoretical predictions for the particular example of 
the "alternating sine flow" and a monochromatic source. We considered a large range 
of values for p, covering more than three orders of magnitude: from p ^ 1 (in which 
case homogenization theory becomes relevant) to p ^ 1 (in which case diffusion starts 
to dominate). The theoretical results were shown in figure [2] to give a qualitatively good 
description of the non-trivial dependence of kd on p. In particular, the numerical results 
were found to share a similar scaling behaviour with the diffusion-dominated Regime I, 
saturating it for large values of p, and the Batchelor Regime III, with the agreement for 
the latter being closer for larger values of Pe. The scaling of Regime IV did not appear 
in our numerical results. Instead, we found that for sufficiently small p, the numerical 
results match prediction (|5.7p that corresponds to Regime V. Thus, for these values 
of p, homogenization theory provides a better estimate to the estimates derived from 
the bounds. Since the transition between Regimes IV and V depends on the value of 
the Strouhal number St (see equation (|5.8p and figure [3]) , we anticipate that when this 
becomes sufficiently small, the range of validity of Regime IV will become sufficiently 
large for its scaling to be realized in the numerics. 

The numerically obtained values were, in some cases, found to be more than one order of 
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magnitude smaller than the values estimated by the bounds (the exception being Regime 
I which for p 3> 1 is saturated). An enhancement of the bo unds can, in general, be mad e 
by finding the opt in ial 'test field' ^ in c onstra int p.ip (see lDoering fc Thiffeault ( 2006 ): 
Shaw et al\ (|2007l) : IPlasting fc YounsJ (|2006l )). Such variational methods are expectec. 
to only im prove the prefactor and not the scaling of bound (|3.4p (and thus Regime 
IV). In fact iPlasting fc Yound (|2006[ ) have shown that for the particular example of the 
"alternating sine flow" and monochromatic source, the choice "0 = s is optimal and any 
improvement relies on knowledge of x that is generally unknown. Conversely, the leading- 
order term in bound p.l3p is independent of the choice of the 'test field' and thus, with 
the current constraints, such methods are not likely to improve the bound in Regime 
III. It should be noted, however, that the biggest advantage of these upper bounds lies 
in predicting (or, to be more exact, restricting) the scaling behaviour of the dissipation 
wavenumber and how this is controlled by the parameters of the system. In that respect, 
the present investigation has proven to be very fruitful. 

We now discuss the relation of our results with a partic ular set of mixing measures , 
the so-called mixing efficiencies, Ep. These were defined in lDoering fc Thiffeault (l2006l) 
and , Shaw et at. (.2007,1 in terms of (1^7^6*1^), for p € Z, and the same variances obtained 
for 6q satisfying (|2.ip in the absence of a flow (u = 0): 



(|VP6l|2) ■ 



(6.1) 



Ep arc commonly larger than unity. In the high-Peclct limit and for spatially smooth 



source flclds, they were shown to satisfy Ep < Pe/p, for p = - 
and equation (|2.9p into the definition for Eg, we obtain that 



Eq = C4 



P/2 



1,0,1. Using 6*0 = ^V-^s 



(6.2) 



where C4 is a non-dimensional number defined as C4 = 

((A-is)2)i/2/5' (the hat symbol 
denotes differentiation with respect to x = x/£s). Similarly, using equation (|2.6|) . 

kd^s 



Ei^C5 



(6.3) 



where Cs is a non-dimensional number defined as cs = 

((V-is)2)V2/5'. Thus, neither 
Eq nor El include separate information about or £,^1- Since we have no control over 
the value of ^g.s, we cannot directly compare the bounds for kd with those for Eq or Ei. 
Instead, the two sets of bounds provide complementary information. We note that from 
equations ()6.2|) and (j6.3|) . wc expect that if the suppression of variance is solely due to 
the suppression of (,e,s (the case of a uniform fiow), the two efficiencies Eq and Ei should 
scale similarly with Pe. If however the suppression of the variance is due to an increase in 
kd, Eq and Ei are expected to scale differently. A separate investigation of the behaviour 
of kd and £g_s will be useful to clarify the types of flow that suppress the scalar variance 
mainly due to transport and those that do so mainly due to stirring. 

Throughout this paper we have been working under the assumption that the spatial 
gradients of the velocity and the source fields are finite. Still, it is worth speculating on 
the implicati ons of our results for rough sources and flows. Th e case of rough sources was 
considered in iDoering fc Thiffeault! (|2006l ) ; IShaw et ali (|2007 ) for the mixing efficiencies 
(|6.ip . In this case, the roughness exponent of the source becomes crucial. For our bound 
(|3.13p , the source roughness will change the balance of the three terms inside the square 
root in equation p.l3p . giving rise to different scalings for the dissipation lengthscale. A 
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detailed examination would need to be performed to detemine the scaling behaviour in 
this case. 

The case of rough velocity fields is also very important because of its relevance to 
turbulent flows. Although in this case bound p.l3|) can not be defined, it is still worth 
examining the implications of our results using simple scaling arguments at the cost of 
losing some of the mathematical rigor. In three-dimensional turbulent fiows, the most 
energetic scales, Lf, are large and control the transport between the sources and the sinks. 
Conversely, the smallest eddies have the largest shearing rate and control the stirring. A 
transition in the behaviour of kd is thus expected when bound p.4p (that is still valid for 
rough velocity fields) intersects the Batchelor scaling, kd ^ i^^. In terms of the Reynolds 
number. Re = ULf /v (where v is the kinematic viscosity), the Batchelor scale reads, 

^'.-^Re-^ (6.4) 

where we assume that Pc Re. Comparing (|6.4p with bound (|3.4p . we obtain that a 
transition occurs when Is ''^ t*s = LfRe~^ . For £s ^ ^s, the scaling k'^ ^ holds while 
for is ^ il, bound (|3.4[) becomes smaller than and the scaling induced from (|4.4p is 
expected. This prediction, however, should still be verified by numerical simulations. 
We plan to address a number of the above mentioned issues in our future work. 
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